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ABSTRACT 

Aims. The K2 mission has recently begun to discover new and diverse planetary systems. In December 2014, Campaign 1 data from 
the mission was released, providing high-precision photometry for -22000 objects over an 80-day timespan. We searched these data 
with the aim of detecting more important new objects. 

Methods. Our search through two separate pipelines led to the independent discovery of K2-19b & c, a two-planet system of Neptune¬ 
sized objects (4.2 and 7.2 R e ), orbiting a K dwarf extremely close to the 3:2 mean motion resonance. The two planets each show 
transits, sometimes simultaneously owing to their proximity to resonance and the alignment of conjunctions. 

Results. We obtained further ground-based photometry of the larger planet with the NITES telescope, demonstrating the presence 
of large transit timing variations (TTVs), and used the observed TTVs to place mass constraints on the transiting objects under 
the hypothesis that the objects are near but not in resonance. We then statistically validated the planets through the PASTIS tool, 
independently of the TTV analysis. 

Key words, planets and satellites: detection, dynamical evolution and stability, individual: K2-19b, individual: K2-19c, general 


1. Introduction 


With the steady release of data from the K2 satellite, several 
projects have begun to search for previously undiscovered plan¬ 
etary systems. A number of interesting systems have already 
come to light ([Crossfield et al.||2015[ |Vanderburg et al.||2014[ 
Foreman-Mackey et al. 2015]). For these systems we now have 
photometry that approaches the precision of the Kepler prime 
mission, and crucially, of host stars much brighter than the typ¬ 
ical Kepler case. This promises the use of radial velocity and 
other techniques to add to our knowledge of these already in¬ 
teresting objects. In this work we present a two-planet system 
observed in K2 field 1. This system, K2-19 (EPIC201505350, 
RA: 11:39:50.476, DEC: +00:36:12.87, Kepmag 12.8), lies ex¬ 
ceptionally close to the 3:2 mean motion resonance (MMR), so it 


* Using observations made with SOPHIE on the 1.93-m telescope at 
Observatoire de Haute-Provence (CNRS), France. 


has the potential to show particularly large TTVs (a concept first 
suggested by | Agol et al.] ( 2005) ; |Hol man ( 2005) for the general 


case). In terms of period ratio, only one object is as yet known 
to be closer to this resonance (and does not show TTVs, due to a 
large libration period). K2-19 was originally presented as a can¬ 
didate planetary system in [Foreman-Mackey et al.| ( [26l5| ) and is 
validated here using further observations. 


The 3:2 MMR is especially significant in both solar system 
and extrasolar planetary systems. For decades, Pluto and Nep¬ 
tune were classified as the only solar system resonant planet pair, 
and their orbits evolve inside of a 3:2 MMR. The Grand Tack 
model, a scenario proposed to explain the current architecture of 
the inner solar system, asserts that Jupiter and Saturn were once 
captured i nto a 3:2 MMR while embedded in t he nascent so¬ 
lar nebula ( |Walsh et al.|2011[|Pierens et al.|2014| ). Furthermore, 
the first extrasolar planetary system ever confirmed, around PSR 
1257+12 ( [Wolszczan & Frail|1992[|Wolszczan|1994| , includes 
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two planets whose orbits are t ightly coupled and are very close to 
residing within the 3:2 MMR ( Malho tra et al.|1992t|Gozdziewski| 


|et al. |2005[ [Callegari et sA.\2W6f 

Transiting exoplanets orbiting main sequence stars represent 
the majority of known planets, but usually lack the necessary 
constraints to allow one to definitively assign membership to 
an individual MMbQ A commonly-used definition of MMR be¬ 
tween two planets is the resulting dynamical state when a partic¬ 
ular linear combination of mean longitudes, longitudes of peri- 
centre, and sometimes longitudes of ascending node librate (or 
oscillate) about 0° or 180° over a given time interval (e.g. 
ray & Dermott| 1 999) . For K2-19 and the 3:2 MMR, this 


investigated (Papaloizou & Szuszkiewicz 2005; Hadjidemetriou 
& Voyatzis||2010| |Emeryanenko||2012 Qgi hara & Kobayashi 


Mur- 


combination is represented by either the angle 6 \ or 62 , where 


0i = 3 A c - 2A b - rn b 
6 2 = 3 A c - 2A b - m c . 


( 1 ) 

( 2 ) 


with A the mean longitude and w the longitude of pericentre. 
Here, and throughout the paper, we label the inner planet as b 
and the outer as c. This widely-adopted definition has several 
shortcomings, which are outlined in Section 1 of Petrovich et 


[ak|p013| . A definition which overcomes these shortcomings is: 
the term resonant is a trajectory that evolves in the region of 
phase space surrounded by the separatricies of a given integrable 
system Hamiltonian (see, e.g. |Morbidelli|2002| ). 

Because time series of planetary mean longitudes and longi¬ 
tudes of pericentre are typically not available in extrasolar sys¬ 
tems, a common practice is to use orbital period ratios by them¬ 
selves as a proxy for resonance. High frequencies of systems 
with ratios of 1.5 and 2.0 (Lissauer et aT7| |201 ip are suggestive 
that the 3:2 and 2:1 commensurabilities represent a significant 
tracer of formation, regardless of whether those planetary can¬ 
didates are actually locked inside of a MMR. Recently, the fre¬ 
quency of systems just outside of the 3:2 period commensura- 
bility has exhibited a distinct asymmetry (e.g. Fabrycky et al 


2014]), which has led to substantial theoretical scrutiny (Batygin 


& Morbidelli 

2013] 

Lee et al. 2013[[Petrovich et al. 2013; Chat- 

terjee & Ford 

2014; 

Delisle et al. 2014;; Delisle & Laskar 2014). 


The period ratio of K2-19b & c as displayed in the K2 data is 
1.5035 14 ^ 0000057 , amon S ^e closest systems to a 3:2 commen- 
surability so far detected. 

We searched the Exoplanet Orbit Database ( |Han et al.|2014[ ) 
for other systems close to this commensurability. The only sys¬ 
tem we could find with a closer normalised distance to resonance 
(defined in|LJt hwick et al. 2012|, A, was the Kepler-372cd pair 
( [Rowe et al.|2014| ), where A is -0.0003, as compared to K2-19 


inear 


with -0.0023. However, neither Kepler-372c nor d exhibit TTVs 
during the Kepler observations due to a particularly long pre¬ 
dicted TTV libration period, -70 years. Also worth noting is the 
Kepler-342cd pair ( [Rowe et al.|2014| ), with a A of -0.0027, which 
also does not show TTVs due to a longer libration period. 

Interest in these special period ratios is motivated by both 
the possibility of making deductions about a system’s formation 
channel and its long-term future stability. Convergent migration 
in protoplanetary discs is an effective and popular MMR forma¬ 
tion mechanism (Snellgrove et al. 2001), an d can also achieve 
three-body resonances ( [Peale & Lee 2002| |Libert & Tsiganisj 
201 1 |, although s ome MMRs are harder to lock into than others 
(Rei n et al.|2012 |Tadeu dos Santos et al.|2015| ). Forming the 3:2 
MMR in particular through energy dissipation has been widely 


1 With limited constraints, one can more easily exclude systems from 
existing within MMRs (Veras & Ford 2012| ). 


2013[|Wang & Ji|2014||Zhang et al.|2014| ). Capture into MMRs 

through gravitational scattering alone - after the dissipation of 
the protoplanetary disc - occurs relatively rarely ([Raymond et| 
|al.|2008| . 

In this work we characterise this multi-planet system, vali¬ 
date its planetary nature using observed TTVs and the PASTIS 
tool, and discuss the implications these have for future observa¬ 
tions. 


2. Observations 

2.1. K2 

Observations were made with the Kepler satellite as part of the 
K2 mission between BJD 2456811.57 and 2456890.33, span¬ 
ning -80 days. The K2 mission (Howell et al. |2014| ) is the sur¬ 
vey now being conducted with the repurposed Kepler space tele¬ 
scope, and became fully operational in June 2014. It is surveying 
a series of fields near the ecliptic, returning continuous high- 
precision data over an 80 day period for each field. Despite the 
reaction wheel losses that ended the Kepler prime mission, K2 
has been estimated to be capable of 80ppm precision for V=12 
stars, close to the sensitivity of the primary mission. All data 
will be public, although at the time of writing only campaigns 
0 and 1 have been released. Targets are provided by the Ecliptic 
Plane Input Catalogue (EPIC) which is hosted at the Mikulski 
Archive for Space Telescopes (MAST^] along with the available 
data products. 

Targets in K2 often display significant pointing drift over the 
K2 observations, typically on a timescale of 6 hours on which 
the spacecraft thrusters are fired. This leads to a major source of 
systematic noise in the lightcurve ( [Vanderburg & Johnson|2014] ), 
the removal of which is the key part of our detrending method. 
The full method is explained in [Armstrong et al.| ( |2015| ), but 
is summarised here for clarity. Initially a fixed ape rture, of ra- 
dius 4 pixels in this case and shape as described in | Armstrong) 
|et al.| ( |2015| ), was centred on the brightest target pixel. A raw 
light curve is extracted directly from this aperture, with back¬ 
ground subtracted using the median out-of-aperture pixels. Row 
and Column centroid variations are found for the time series. At 
this point points associated with spacecraft thruster firings are 
removed, and the remaining points decorrelated from the cen¬ 
troid variations. We decorrelate from both row and column cen¬ 
troids simultaneously, as they are not statistically independent. 
This process can leave systematic or instrumental noise in place 
(see For eman- Mackey et al .|2015| ). In this case this noise seems 
to be weak compared to the intrinsic stellar variability, which 
occurs with a magnitude of -1%. 

This stellar variability is removed, along with any longer pe¬ 
riod systematics, through the application of an iteratively fitted 
polynomial. A 3D polynomial is fit to successive 2 day wide 
regions, with the fit repeated for 20 iterations clipping points 
greater than 3 cr from the best fit line at each iteration. This fit is 
then used to detrend a 5 hour region at its centre, and the process 
repeated for each 5 hour region. As the principal components of 
the instrumental noise show variations on order 10+ days in cam¬ 
paign 1 (Foreman-Mackey et_al. 2015[ ) they should be removed 
through this process, and crucially will not affect the transits. We 
found that for this target the best results were obtained by per¬ 
forming this polynomial flattening immediately after extracting 


https ://archive. stsci.edu/k2/ 
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2.3. SOPHIE 



Fig. 1. Data taken from the NITES telescope. The best fit transit derived 
from K2 observations is also plotted, along with the binned light curve, 
and fits the shape seen with NITES well. 


We observed the star K2-19 with the SOPHIE spectrograph 
mounted on the 1.93 m telescope at the Haute-Provence Obser¬ 
vatory (France). We used the high-efficiency mode which has a 
spectral resolution of about 39 000 at 550nm. This mode is pre¬ 


ferred for the observations of relatively faint stars (e.g. Santeme 
|et al.||2014| ). For more information about the SOPHIE spectro¬ 
graph, we refer the interested readers to |Perruchot et al. ([2008 ) 
and Bouchy et al.| ( |2009] ). We secured five epochs between 2015- 
01-23 and 2015-02-02 as part of our on-going TRANSIT consor- 
tiu nf] The exposure time ranges between 400s and 2700s, which 
lead to signal-to-noise ratio per pixel in the continuum at 550nm 
ranging between 13.5 and 22.1. 

We cross correlated the spectra with a numerical mask cor¬ 
responding of a G2 dwarf (Bara nne et al.|1996}|Pepe et al.|2 002) 
and find a unique line profile with a width compatible with the 
rotational period found in K2 photometry (see next section). The 
derived radial velocity has a median uncertainty of 16ms -1 , and 
shows no significant variation. 


the lightcurve, before decorrelating the flux from the centroid 
motion. The resulting light curve is shown in Figure[2] Note that 
some transits of each planet occur simultaneously with the other 
through the dataset. The final such simultaneous transit is of in¬ 
creased duration, and its shape begins to show signs of both plan¬ 
ets individually. 

2.2. NITES 

Due to its proximity to MMR, K2-19 had the potential to show 
significant TTVs (see Section [5]). As such we scheduled it for 
further ground based transit observations. The Near Infra-red 
Transiting ExoplanetS (NITES) Telescope is a semi-robotic 0.4- 
m (f/10) Meade LX200GPS Schmidt-Cassegrain telescope in¬ 
stalled at the ORM, La Palma. The telescope is mounted with 
a Finger Lakes Instrumentation Proline 4710 camera, contain¬ 
ing a 1024 x 1024 pixels deep-depleted CCD made by e2v. The 
telescope has a FOV and pixel scale of 11 x 11 arcmin squared 
and 0.66" pixel -1 , respectively and a peak QE> 90% at 800 nm. 
For more details on the NITES Telescope we refer the reader to 
|McCormac et al.| ( |2014| ). 

One transit of K2-19b was observed on 2015 Feb 28. The 
telescope was defocused slightly to 3.3" FWHM and 814 images 
of 20 s exposure time were obtained with 5 s dead time between 
each. Observations were obtained without a filter. The data were 
bias subtracted and flat field corrected using PyRAI0 and the 
standard routines in IRAFq and aperture photometry was per¬ 
formed using DAOPHOT (Stetso n) 1987| ). Ten nearby compari¬ 
son stars were used and an aperture radius of 6.6" was chosen 
as it returned the minimum root mean square (RMS) scatter in 
the out of transit data. Initial photometric error estimates were 
calculated using the electron noise from the target and the sky 
and the read noise within the aperture. The data were normalised 
with a first order polynomial fitted to the out of transit data. The 
resulting lightcurve is shown in Figure [T] 

3 PyRAF is a product of the Space Telescope Science Institute, which 
is operated by AURA for NASA. 

4 IRAF is distributed by the National Optical Astronomy Observato¬ 
ries, which are operated by the Association of Universities for Research 
in Astronomy, Inc., under cooperative agreement with the National Sci¬ 
ence Foundation. 


2.4. AstraLux 


In order to search for neighbouring stars which could be pro¬ 
ducing the transit signals, we obtain a high-spatial resolution 
image of K2-19 by applying the lucky-imaging technique with 
the AstraLux instrument installed at the 2.2m telescope at the 
Calar Alto Observatory (Almeria, Spain). We obtained 65000 
frames with an individual exposure time of 80 milliseconds in 
the SDSSz band. The images were reduced by using the obser¬ 
vatory pipeline (see |Hormuth et al.|2008| ), which performs a ba¬ 
sic reduction and selects the 10% of frames showing the largest 
Strehl ratios ( [Strehl|p~902] ). These frames are then aligned and 
stacked to produce the final high-spatial resolution image. The 
sensitivity curve of this image was obtained by adding artifi¬ 
cial sources with different magnitude contrasts at different angu¬ 
lar^ separations and measuring the recovery rate (see [Lillo-Box| 
et al.|2014] for further details). We reach contrast magnitudes of 
A m z > = 5.5 mag at 1" and A m Z ' = 6.3 mag for angular sepa¬ 
rations larger than 1.5" at the 5 cr level. No sources were found 
closer than 6" within our sensitivity limits. The analysis of the 
Kepler centroids also revealed no significant shifts during tran¬ 
sits above the Kepler pixel size (~ 4"). 


3. Stellar parameters 

We obtained the parameters of the host star from the spectral 
analysis of five co-added SOPHIE spectra. First, we subtracted 
from the spectra pointing to the source (in fibre A), any sky con¬ 
tamination using the spectra of fibre B, after correcting for the 
relative efficiency of the two fibres. The final spectrum has a S/N 
of the order of 25 around 6070A. 

To derive the atmospheric parameters, namely the effective 
temperature (r ef f), surface gravity (logg), metallicity ([Fe/H]), 
and microturbulence (v m i c ), we followed the methodology de¬ 
scribed in |Tsantaki et al.|p013| . This method relies on the mea¬ 
surement of the equivalent widths (EWs) of Fei and Fen lines 
and by imposing excitation and ionisation equilibrium. The anal¬ 
ysis was performed in local thermodynamic equilibrium (LTE) 
using a grid of (Kuruc z~fl993[) m odel atmospheres and the radia¬ 
tive transfer code MOOG ( |Sneden|1973| ). Due to the low S/N of 

5 OHP program: 14B.PNP.HEBR 
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Table 1 . Stellar parameters for K2-19 


Parameter 

Value 

Units 

7’cir 

log g * 

Unic 

[Fe/H] 

^rot 

5230 + 417 
4.39 ± 0.79 
0.92 ± 0.5 
0.38 ± 0.23 
20 . 3 ^ 

K 

dex 

kms -1 

dex 

days 

Derived Parameters: 



R+ (spectroscopic) 

1.03 + 0.2 

Re 

M* (spectroscopic) 

0.92 + 0.14 

M 0 

log g* (transit) 

4.52 + 0.22 

dex 

R+ (transit) 

0.88 + 0.06 

Re 

M* (transit) 

0.89 + 0.06 

M 0 


our spectrum, the EWs were derived manually using the IRAF 
splot task. 

From the above analysis, we conclude that the host star is a 
slightly metal-rich K dwarf. The derived parameters are shown 
in TableQ] The stellar radius and mass were derived from the cal- 
ibration of|Torres et al.| ( |20 10a), updated with the version from 
Santos et al. ( |2013| ) and the atmospheric parameters described 
above. We also included in Table Q] the determination of surface 
gravity from the transit fit parameters (see Section]?]) and the re¬ 
spective results of stellar mass and radius. We proceed using the 
transit derived parameters, as they are much more accurate. 

We also study the stellar variability inherent in the lightcurve 
of K2-19. This lightcurve may be contaminated by remnant in¬ 
strumental noise, but we find that repeating patterns apparent 
across the entirety of the K2 observations do not generally match 
the principal noise components seen (see |Foreman-Mackey et al.| 
2015, for these components). A weighted, floating mean Lomb- 
Scargle (FS) periodogra m (|Lomb||1976[ Scargle||1982| ), follow¬ 
ing the method of Press & Rybicki (1989), gives a principal pe¬ 
riod of 20.3 days (with an FS statistic of 323000, significantly 
above the background). Assuming this peak is due to stellar ro¬ 
tation then we are able to derive P rot as shown in Table [T] Errors 
on P rot are derived from the FWHM of the periodogram peak. 


4. Light curve fitting 

To obtain the transit shape and parameters we limit ourselves 
to the K2 data, as it is of significantly higher precision and the 
NITES data do not sho w the full transit. The data were detrended 
as described in Section [2T] then cut so that only data within a 7 
transit width region centred on each transit were used. We also 
removed all simultaneous transits, along with two specific points 
in separate transits which showed clear evidence for being within 
a spot crossing (significant brightenings within transit relative to 
their local transit shape). These points are highlighted in Figure 
[2] We note that there are other apparently bad points which were 
not removed - the decision to remove a point was based entirely 
on whether it was clearly anomalous within its local transit, to 
avoid excessive bias, and so points which only appear bad when 
shown phase folded and against the fit will remain. The data were 
then fit using th e JKTEBOP code (e.g. |Southworth|20l3}|Popper| 
|& Et zel 1981), with numerical integration used to account for 
the long cadence of K2 observations (splitting each point into 60 
integrated sub-points). 

We initialised the fits with a linear limb darkening coefficient 
of 0.56, suitable for a K dwarf, which was then allowed to vary. 
We then tested for eccentricity, but found no constraint for either 
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object. As such for the remaining tests the eccentricity of both 
planets was set to zero. To derive robust errors we used a Monte 
Carlo process whereby Gaussian observational errors are added 
to each data point and the fit repeated 1000 times, producing 
a distribution of best fits. The medians and 68.27% confidence 
limits are then taken to produce values and errors. 

This process does not account for systematic errors in the 
light curve, such as starspots or correlated instrumental noise. 
These are of particular concern for K2-19; as has been noted 
there is evidence within some transits for spot crossings. In the 
past such crossings have proven useful in modelling starspots 
(e.g. |Barros et al.||2013[ |Beky et ak||2014| ) but here they form 
a source of contamination to our fits. We test for the effect of 
these spots by adopting a prayer bead style residual permutation 
test. In this process, a best fit is acquired, and then the residuals 
of the data to this fit are ‘rolled’ through the dataset, and a fur¬ 
ther best fit acquired each time. Due to the low cadence of K2 
observations, there are not enough points near transit to get a dis¬ 
tribution of parameter values through this method (270 and 183 
tests respectively for planets b & c). However, the prayer-bead 
generated distribution at least allows us to obtain an estimate of 
the systematic effect on our transit parameters. In all cases these 
systematic errors were comparable to or smaller than the Monte 
Carlo generated errors. As such we present final values and er¬ 
rors from the Monte Carlo tests. While we acknowledge that this 
method of testing for systematic errors is merely an estimate (es¬ 
pecially as the full effect of spots only appears in transit), we 
note that as the errors generated by the prayer-bead process are 
not significantly larger that those from the Monte Carlo tests, the 
effect of systematics on the transit parameters is not particularly 
strong. 

The resulting best fits are shown for each planet in Figure 
[2] Note that the derived ephemeris are taken from only a small 
part of the TTV phase curve and so will require correction; see 
Section [7] for detail. In particular, there are significantly larger 
errors on the period when TTVs are taken into account - final 
values are found in Section 17^21 


5. Transit timing variations 

Given the periods found in Section [4] it is immediately appar¬ 
ent that the two planets in the K2-19 system lie close to the 3:2 
MMR. It is common for systems close to MMR to show partic¬ 
ularly large TTVs ( [Lithwick et al.||2012[ |Xie||2014| ) and so we 
searched the data in the hope of seeing variations. This search 
was carried out using the transit shape defined by our best fit pa¬ 
rameters. As before, simultaneous transits were ignored. Points 
marked previously as being clear spot crossings were also ex¬ 
cluded. We cut the data to a region within 2 full transit widths 
of the approximate transit centres, then passed the model transit 
over this region with a resolution of 0.00015 days. The mini¬ 
mum x 1 of this test series was recorded, at which point each 
datapoint was perturbed by a random Gaussian with standard 
deviation equal to the point error. The fit was then repeated, and 
this process undergone for 1000 iterations, to get a distribution 
of transit times. The mean of the distribution is then taken as the 
transit time. As when fitting the transit shape, this process does 
not account for systematic errors. This is particularly concerning 
for measuring transit times, because due to the low cadence of 
K2 observations only a few points are seen within each transit. If 
one of these points is significantly perturbed by a spot crossing 
(which occurs visibly for some transits) then the measured time 
would be strongly affected. To estimate the effect of these sys¬ 
tematics, we repeat the prayer-bead residual analysis of Section 
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Phase Phase 


Fig. 2. Top: Extracted light curve for K2-19, showing significant stellar variability. Middle: Flattened and detrended lightcurve, showing transits of 
the inner planet (b, red), outer planet (c, green) and simultaneous transits of both planets (magenta). Some outlier points are not shown for clarity. 
Bottom Left: The phase folded transits of planet b, excluding simultaneous transits and showing the best fit model. Some points (shown lighter 
than the others) displayed clear evidence of spot crossings by the planet and were excluded from the fit. Bottom Right: Same for planet c. Note the 
change in y-axis scale. 

Table 2. System parameters 


Parameter 

Units 

b 

c 

Model Parameters: 

P 

days 

7 919454+ 0 - 000081 

/.yiy^j tt _ 0>000078 

11 90701+ 000039 

ii.^u/ui _ 000044 

To 

BJDj-db - 2456000 

qiq o qq /1 c+0.00036 
o 13.3 o34 J _ 0 00039 

817.2759 + 0.0012 

Rp/R* 


0 0753 +0 - 0028 
U,U/JJ _ 0>0015 

0 0439 +0 0011 

w.wh-o ^_ 00012 

( R p + R*)/a 


0 0572 +0 - 0084 

U,UJ/ - 0 . 0042 

0 0414+ 0 - 0015 
w.wh-ih -_ 00009 

i 

deg 

QO oq+1.08 

oo.o3 _ 0 89 

89 91 + 0.05 
oy ’ yL -0.32 

e 


0 (adopted) 

0 (adopted) 

Limb-Darkening 


0.552 ± 0.041 

0 57 +0 - 14 

Derived Parameters: 

R P 

^© 

7 oq+0.56 
-0.51 

4.21 ±0.31 

a 

AU 

0.077” 

0 1032 +a0074 
U,1UJ - 0 . 0080 

R inc 

S © 

87-7!” 

48.8!^ 

Pc/Pb 


1 503514 +0,0000 '^ 

i.JUJJl‘t _ 0 000057 


A 


0.00234 ± 0.00002 



Notes. A is defined in Section]?] and represents the normalised distance to resonance. Note that Pfc,P c , and parameters derived from them are only 
instantaneous measurements, and will change over the course of the TTV phase curve (see Section [TJ. Transit based stellar parameters are used 
for derived quantities. 


[4] In this case though, as we are considering each transit inde¬ 
pendently, there are even fewer data points near transit (typically 
~30). Also of concern is that the full effect of spots can gener¬ 
ally only be seen when they are occulted in transit, where there 
are even fewer points. As such we perform this analysis and es¬ 
timate the systematic contribution to our error budget by taking 
the maximum and minimum parameter values which arise from 


the prayer bead test, over the ~30 iterations. We adopt these most 
pessimistic values as our lcr errors, to ensure that we do not un¬ 
derestimate the errors on our transit times. The adopted values 
(from the mean of the monte carlo distribution) and errors (from 
the maximum and minimum of the prayer-bead residual test) are 
given in Table [3] 
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Table 3. Detected transit times 


Planet Time (BJD tdb - 2456000) Error Source 


b 

813.3841 

0.0016 

K2 

20 

b 

821.3039 

0.0107 

K2 


b 

837.1382 

0.0014 

K2 

15 

b 

845.06176 

0.00098 

K2 


b 

860.9000 

0.0012 

K2 

B !0 

b 

868.8196 

0.0016 

K2 

13 

C 

b 

884.6597 

0.0017 

K2 

E 

U 5 

b 

1082.6895 

0.0022 

NITES 

6 


C 

817.2741 

0.0032 

K2 

c 

841.0942 

0.0068 

K2 

c 

864.9105 

0.0069 

K2 

c 

888.7136 

0.0038 

K2 


Notes. Simultaneous transits are not shown here. 


3 4 5 6 

Orbit 



Fig. 3. Observed-Calculated (O-C) transit times for planet b. Calculated 
times are taken by setting the O-C values at the first and last transits 
to zero because of the absence of well-determined period information. 
We do not use the ^-derived ephemeris as these may misrepresent the 
scale of the TTYs. 


We were fortunate enough to obtain an additional transit of 
planet b with the NITES telescope (Section |Z2| ). The time of 
this transit was obtained using the transit shape derived from the 
K2 observations, meaning the only fit parameter was the time 
of transit centre. The same Monte Carlo test was performed, but 
now adopting the standard deviation of the resulting distribution 
as the error. We did not repeat the prayer-bead analysis in this 
case, as the transit did not show evidence for systematics. 

The observed-calculated times found from our TTV analy¬ 
sis are shown in Figures [3] and [4} calculated as described in the 
Figure captions. Planet b, in particular the NITES observation, 
show large TTVs of over an hour from the expected K2, linear 
ephemeris based, time (~a quarter of the transit duration). Within 
the K2 data alone we do not find the TTVs to be significant, be¬ 
yond the third non-simultaneous transit for planet b which ar¬ 
rives earlier than would be expected. An initial analysis of these 
TTVs is performed in Section |T2| We note that this detection 
of TTVs implies that the ephemeris in Table ^are likely not the 
‘true ’ ephemeris , in the sense of the mean transit interval over 
long timeframes. Readers should thus be careful in predicting 
transit times. This is discussed further in Section [ 7 ] 


Fig. 4. As Figure [ 5 ] for planet c. In this case the calculated times are 
taken from the ephemeris of Table [2] 


6. PASTIS validation 

6.1. Overview 

Candidates from Kepler or K2 can be confirmed or validated in 
a number of ways. These include radial velocity observations, 
BLENDER and PASTIS analyses of t he probability of false posi¬ 
tive scenarios ( [Torres et al.||2010b| |Diaz et aT7||2Q14| ), high res¬ 
olution photometry to search for close companions (e.g. [Everett 
|et al.|2015|[Lillo-Box et al.|2012[|2Q14[|Law et al.|2014| ) or TTVs 
(e.g. [Steffen et al.|2012| Nesvorny et al.|2014|). 

Astrophysical false-positive scenarios such as eclipsing 
binaries might mimic the transit of a planet (Torres eTal. j2005| . 
Before claiming the planetary nature of a small and periodic 
signal, one should first rule out the possibility that this signal 
has a non-planetary origin. When two or more sets of transits 
are detected on the same target, their probability not to be 
planets significantly decreases ( Lissauer et al.||20T2} |Lissauer 
|et al.||20T4] ). Rowe et al. ( |2014| ) used this planet-likelihood 
“multiplicity boost” to validate a large sample of planets in 
multiple systems. This validation was possible because of the 
large sample of systems considered: a significant fraction of 
these systems are expected to be planets. 

This argument cannot be used to validate a single multi¬ 
planetary system. Even if the two sets of transit in the K2-19 
system are close to MMR, they could orbit two different stars. 
The period commensurability between the two could be the re¬ 
sult of chance. For example, [Lissauer et al. ( |2014| ) presented the 
interesting case of Kepler-132, a binary system in which each 
stellar member hosts a transiting planet at 6.2 and 6.4 days, re¬ 
spectively. That these two planets are orbiting with almost the 
same period is likely a coincidence. Therefore, the fact that the 
target star K2-19 hosts two planet-like objects close to MMR 
might also be a coincidence, so it is not evident that they are 
planets and that they both transit the same star, even if this is 
quite likely. In this section, we verify whether both transit sets 
have a planetary origin and orbit the target star. 


6.2. Definition of the scenarios and general framework 

We identified 11 scenarios that could explain the K2 light curve, 
which we called So to S\o’ 
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- So\ the target star K2-19 is transited by two planets; 

- S i: a planet is transiting the target star K2-19 every 8 days 
and another planet is transiting a physical companion to the 
target star every 12 days; 

- S 2 '. a planet is transiting the target star K2-19 every 12 days 
and another planet is transiting a physical companion to the 
target star every 8 days; 

- £ 3 : two planets are transiting a physical companion to the 
target star; 

- S 4 : a planet is transiting the target star K2-19 every 8 days 
and an eclipsing binary with a period of 12 days is bound 
with the target star; 

- S 5 '. a planet is transiting the target star K2-19 every 12 days 
and an eclipsing binary with a period of 8 days is bound with 
the target star; 

- So', a planet is transiting the target star K2-19 every 8 days 
and an eclipsing binary with a period of 12 days is chance- 
aligned with the target star (background / foreground); 

- S 7 '. a planet is transiting the target star K2-19 every 12 days 
and an eclipsing binary with a period of 8 days is chance- 
aligned with the target star (background / foreground); 

- S&: a planet is transiting the target star K2-19 every 8 days 
and another planet is transiting every 12 days a star chance- 
aligned with the target star (background / foreground); 

- Sg\ a. planet is transiting the target star K2-19 every 12 days 
and another planet is transiting every 8 days a star chance- 
aligned with the target star (background / foreground); 

- S 10 : two planets are transiting a star chance-aligned with the 
target star (background / foreground). 

To he lp th e reader, a sketch of these scenarios is displayed 
in Figure |A.1| The scenario So is the one we want to test. The 
scenarios S 1 to S 5 invoke a stellar companion to the target star 
which is orbited either by planets or by eclipsing binaries. The 
scenarios S 6 to <Sio invoke a chance-aligned star transited by 
planets or eclipsed by lower-mass stars. For stability reasons, 
we did not consider here the case where two stars are eclipsing 
the same star at 8 and 12 days. 


The probability of a scenario Si is defined in the Bayesian 
framework as P ( Si | D, J), with P the probability, D the data, 
and I the information. Using Bayes theorem, this probability 
can be expressed as: 


P(Si\D,D 


P(Si\DxP(D\Sj,D 

P(D\I) 


(3) 


P (Si 1 1) is the a priori probability of the scenario S» 
P(D\Si,I) is its marginalised posterior probability, and 
P (D1 1) is a normalisation factor. The latter one being difficult 
to estimate, we instead computed the odds ratio Oij between two 
scenarios S t and Sj, so that this normalisation factor is cancelled 
out: 

c _ P(Sj\D,D _ P(Si\DxP(D\Si,D 

IJ p(Sj\v,1) r(Sj\i)y.r(v\Sj,i)' 

Then, assuming that ^=0 ^ I -0 = 1, the probability of a 
given scenario Si can be computed as: 


P(Si\D,D = 


( 10 
yj =0 



(5) 


Therefore, to estimate the probability of the scenario So , we need 
to compute the odds ratio between each pair of scenarios. For 


that, it is necessary t o ev aluate the a priori probability of each 
scenario (see Section [63]) as well as their marginalised posterior 
probability (see Section |6.4| ), which is defined as: 

P(D\ Si,I) = fp(0\ Si, I)xP(D \e, Si, I) dO, (6) 

with 6 the parameter space needed to model the scenario Si. 
Note that P(D\0, Si , I) is also called the data likelihood which 
we computed assuming the data points are independent and 
normally distributed around their mean value {x\, ..., x n ) with a 
width of {<Ti,..., o~ n }, with n the number of data points: 


f P(D\e,Si,I) = fl 2_ exp 
1 1 V2 ntr k 


k= 1 


1 / x k - m k (6 1 Sj) 

2 \ cr k 


(7) 


with m k ( 6 1 Si), the value of the model describing the scenario Si 
that corresponds to the data point x k . 


6.3. A priori probabilities 

To compute the probability of each scenario using the equations 
[3] to [7J we first need to determine their a priori probability. In 
the literature, there is no estimation of what is the occurrence 
rate of having, e.g. a two-planet system (like <So) or a planet in a 
triple system (like S 4 and S 5 ). Thus, we can not rely on robust 
statistical analysis of these scenarios to use as priors for our 
validation. To estimate these probabilities, we assume that each 
element of the systems has an independent probability to occur, 
i.e. the probability of having a planet is independent from that of 
having a stellar companion. This assumption is obviously wrong 
since stellar multiplicity should affect the formation of planets 
(e.g. | Wang et ah||2014| ), unless the stars are far enough from 
each other so that the impact is small. However, we assume this 
discrepancy does not significantly change our results. 


For the stellar multiplicity, we used the results from 
|Raghavan et al.| ( |2010[ ) who reported an occurrence of binary 
star systems of 34% and an occurrence of triple systems of 
9%. Then, we used the results from Mayor _et al. ( |2011| ) who 
reported an occurrence of planets at the level of 75% (for any 
planets with periods up to 10 years). Assuming independence 
between the occurrence rates, we estimated that the a priori 
probability to have, e.g. a target orbited by two planets is 
P(So 1 1) = 0.75 x 0.75 = 0.5625, or the probability to have 
two planets orbiting two different stars of a binary system is 
P (Si | J) = P (S 2 1 1) = 0.75 x 0.75 x 0.34 = 0.19125. We 
listed the a priori probabilities in the Table]?] 

To estimate the a priori probability to have a background 
source of false positive (scenarios S6 to <Sio), we use the As- 
traLux data (Section \2A\ . This data confirmed that no nearby 
star hosts one of the transiting objects. However, the data does 
not completely rule out the possibility that a background or fore¬ 
ground star, chance-aligned with the target, hosts one of the tran¬ 
sits. To evaluate this probability, we used the Besan 9 on galactic 
model ( [Robin et al.|2003| ) and simulated a stellar catalog of one 
square degree around the target star. We considered all types of 
stars with an apparent magnitude between 10 and 20 in the z' 
band. We set the interstellar extinction to zero magnitude per 
kpc in the simulation, which we corrected a posteriori using the 
galactic extinction model of Amores & Lepine ( 2005). We de¬ 
rived the expected density of stars in the vicinity of K2-19 per 
bin of magnitudes within different angular separations from the 
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Table 4. A prio ri probability for the various scenarios given at the start 
of Section lOl 


Probability 

Equation 

Value 

P(SoU) 

0.75 x 0.75 

0.5625 

PCS! ID 

0.75 x 0.34 x 0.75 

0.19125 

P(S 2 \D 

0.75 x 0.34 x 0.75 

0.19125 

rcSiiD 

0.34 x 0.75 x 0.75 

0.19125 

PCSa\D 

0.75 x 0.09 

0.0675 

pcs 5 \d 

0.75 x 0.09 

0.0675 

pcs 6 \d 

0.75 x (4.6 x 10“ 5 ) x 0.34 

1.17 x 10“ 5 

PCSi\D 

0.75 x (2.6 x 10“ 5 ) x 0.34 

6.57 x 10“ 6 

P(Sz\D 

0.75 x (4.6 x 10- 5 ) x 0.75 

2.58 x 10“ 5 

PCS 9 \I) 

0.75 x (2.6 x 10“ 5 ) x 0.75 

1.45 x 10“ 5 

P(S W \D 

(2.6 x lO' 5 ) x 0.75 x 0.75 

1.45 x 10“ 5 


target star. This density of stars is plotted in Figure [5] together 
with the 5-cr detection limits from AstraLux lucky imaging. 

A background source of false positives can only mimic the 
transit depths of the planets if its magnitude is within a range Am 
defined as ( [Morton & Johnson|2011 ): 

Am = 2.51og 10 ||^j, (8) 

where S tr and 6 b g are the depths of the transit, as measured in the 
light curve, and the depth on the background star, respectively. 
Assuming a maximum eclipse depth of 50% for the background 
star, we find that the maximum magnitude range is of 4.9 and 
6.0 for the transit signal at 8 and 12 days (respectively). The 
target star is of magnitude z' = 12.6, hence false positives can 
probe stars as faint as magnitude 17.4 and 18.6 (respectively) in 
the same bandpass. These maximum magnitude limits are also 
represented in the Fig. [5] By integrating the expected number of 
stars that are bright enough to mimic the transit depths if they 
were eclipsing binaries within the detection limits of AstraLux, 
we find 2.6 x 10 -5 stars that might mimic the 8-d signal and 
4.6 x 10 -5 stars that might mimic the 12-d signal. To derive the a 
priori probability of the scenarios So to <Sio, we therefore multi¬ 
ply these values with the probability of being an eclipsing binary 
or that of hosting a planet. These values are reported in Table [4] 
Note that the a priori probability for scenarios to <Sio, where 
the chance-aligned star is transited by one or both planets is over¬ 
estimated since no transiting planet has been found so far with a 
transit depth as large as 50%. 


6.4. The posterior distribution 

To estimate the posterior probability that the K2-19 data are 
produced by one of the aforementioned scenarios, we used 
the PASTIS software ( [Diaz et ak||2014[ |Santeme et al.||2015] ) 


The light curve 


to model the K2 photometric measurement^ 
was modelled using the EB0P code ( [Nelson & Davis|[l972 
|Etzel||198H|Popper & Etzel||1981 ) extracted from the JKTEB0P 


package (Southworth 2008). For the limb darkening coefficient, 
we used the interpolated values fro m|Claret & Bloemen| ( |20lT] ). 
To model the stars (see |Diaz et al.||2014[ for a more detailed 
description on how we model the stars in PASTIS), we used 
the Dartmouth stellar evolution tracks of potter et al.| (|2008| 
and the BT-SETTL stellar atmosphere models of Allard et al. 


6 We used only the data collected in the vicinity of non-simultaneous 
transits. 
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Fig. 5. Map of the density of background stars chance-aligned with K2- 
19, integrated within an angular separation of up to 4", as a function of 
its magnitude in the z'-band. The negative-slope hatched region displays 
all the stars that would have been significantly detected in the AstraLux 
lucky-imaging data with more than 5-cr. The negative-slope and verti¬ 
cal hatched regions display the stars that are too faint to reproduce the 
observed transit depth of K2-19b & c, respectively. 


(2012| which we integrated in the Kepler bandpass. We used the 
results from the spectroscopic analysis for the parameters of the 
target star and the initial mass function from |Kroupa| ( |200 1) for 
the blended stars. The stars that are defined as gravitationally 
bound are assumed to have the same metallicity and the same 
age. The orbits are assumed to be circular. We assumed that the 
blended stars are at least fainter by 1 magnitude in the V-band 
than the target star, otherwise any such star would have been 
clearly identified in the spectroscopic data. The K2 light curve 
was modelled allowing the out-of-transit flux, contamination 
and an extra source of white noise (jitter) to vary. We used an 
oversampling factor of 10, to account for the finite integration 
time of the Kepler long-cadence data ( |Kipping||2010| ). We also 
model the spectral energy distribution of K2-19 composed by 
the magnitudes in the Johnson B and V, Sloan g' and i', 2-MASS 
J, H, and Ks, and WISE W1 t o W3 bandpasses from the APASS 
database ([Henden et al.|2015|) and the A11WISE catalog (|Wright| 
|et al.|2010| ). 


We analyse the aforementioned data using the MCMC 

( |2014| . For the orbital 


procedure described in Diaz et al. 


ephemeris, we used Normal priors matching the ephemeris 
reported in Table [2] with uncertainties boosted by 100, to 
avoid biasing the results with too narrow priors. For the other 
parameters, we choose uninformative priors. We limit the priors 
on the planet radius to be less than 2.2 R j up , which is the 


radius of the biggest planet found to date: KOI-13 (Szabo et 


|al.|2011| . The exhaustive list of parameters and their priors are 
reported in Table [A] For all scenarios, we ran a minimum of 
10 chains of l.Kriterations, randomly started from the joint 
prior distribution. We then selected the chains that converged 
toward the maximum likelihood, which we thinned, to derive 
the posterior distributions for the 11 scenarios. If the posterior 
distributions had less than 1000 independent samples, we ran 
new chains until reaching this threshold. For the scenarios So 
to <Sio, we needed to run up to 40 chains to reach the threshold. 
We report in Table [A] the median and 68.3% confidence interval 
for the free parameters. All the fitted parameters for the scenario 
So are compatible within 1-cr with those derived in Section [4] 
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6.5. Scenario probability and planet validation 

We marginalised the posterior distribution P(6\Si,I) x 
P(D\ 6, Si, I) over the parameter space using the method de¬ 
scribed in |Tuomi & Jones | ( |2012| ). As discussed by the authors, 
this method underestimates the Occam’s razor when computing 
the evidence. However, as already discussed in |Diaz et ah] 
( 2014| ) and |Santerne et al] ( |2014| ), this limitation is expected 
to be relatively weak here since most scenarios have the same 
number of free parameters and most parameters share the same 
priors (e.g. target parameters, orbital parameters). Therefore, we 
assumed our results are not dependent on the method used to 
marginalise the posterior distribution. 


We then computed the probability of each scenario using 
Equation [5] which we multiplied by the a priori probabilities 
listed in Table [4] These scenarios’ probabilities are reported in 
Table [5] The scenario So (i.e. two planets transiting the target 
star) has a probability of 99.2% while the other scenarios have 
a probability below 1%. We can therefore conclude that (1) the 
two objects producing the transits seen in the K2-19 light curves 
are planets, (2) the two planets are transiting the same star, and 
(3) that star is the target star K2-19. As such we are able to sta¬ 
tistically validate the two planets K2-19b & c. 


Table 5. Probability for the various scenarios. Note that the probabili¬ 
ties of the scenarios S\ to <Sio are so low that we provide instead their 
logarithmic values. The logarithmic value of scenario So is provided for 
comparison. 


Probability 

Value 

P(S 0 \D,I) 

99.2!“% 

\og t0 (F(S»\D,I)) 

-0 0035 +00017 

U.UVJJ^_ 00035 

log w CP (SOD,!)) 

-2.17 + 0.30 

l°g 10 (P (<S 2 1 D, J)) 

-2.90 ± 0.32 

logi 0 (IP (S 3 | £>,/)) 

-3.71+0.34 

log lo CP(<S 4 |0,J)) 

-3.80 + 0.31 

log 10 (^(5 5 |D,J)) 

-4.42 + 0.44 

logioCPGSel £>,/)) 

-5.03 + 0.54 

log 10 CP(<S 7 |2>,J)) 

-5.17 + 0.38 

i°g 10 rp (*s 8 1 £>, i)) 

-4.41 + 0.41 

] Og|0 (V(S,\'DJ)) 

-3.87 + 0.39 

l°g 10 CP (*S 10 1D, J)) 

-6.00 + 0.32 


7. TTV mass limits and stability 

7.1. Transit timing variations - overview 

We leave a detailed study of the TTVs to future analysis (Barros 
et al in prep). It is however possible to place a number of con¬ 
straints on the system even with the limited coverage of the TTV 
phase curve which we obtain here. For this initial analysis, we 
use the analytical representation of the TTVs derived by |Lith-| 
| wick et aT7| ( |2() 12 ), hereafter L12, which has been shown to be 
valid for systems near MMR ( [Deck & A gol 2014), a condition 
strongly met in this case. We note that the L12 formulae are only 
valid for objects not explicitly ‘in’ resonance. Without further 
constraints it is impossible to tell whether this condition is met 
here. As such this analysis proceeds under the assumption that 
the objects are not in resonance. The use of the L12 formulae 
allows us to obtain a more intuitive description of the parame¬ 
ter space than is generally possible using N-body simulations. 


Given the potential for spots or other systematic errors to affect 
the K2 transit times, and the otherwise limited coverage of the 
TTV phase curve, we defer such an analysis to future work. 

The TTV phase curve described by L12 is a sinusoid with 
two key parameters: (1) an amplitude |V| given as a function of 
planetary mass, stellar mass, A (the normalised distance to reso¬ 
nance), and the free eccentricity Zf ree (a complex number), and 
(2) a period given by 

p - P ° uter (Q l 

r super — .i | \ y ) 

j|A| 

where 

A = ^out *iz±_ l (10) 

dinner j 

For the 3:2 MMR j - 3. This leads to a phase curve of the 
form 


TTV = \V\sin( t —^+ct >) (11) 

^super 

where 0 is the phase of the curve and changes over the secular 
timescale (hence is constant for our purposes). In our case, we 
can set to to be the time of first transit to acceptable accuracy, 
due to the alignment of planetary conjunctions demonstrated by 
the simultaneous transits observed. Both the amplitude and pe¬ 
riod depend strongly on A. The closer a system is to resonance, 
the larger the amplitude becomes, but the longer the period. For 
a system as close to resonance as K2-19, the period is particu¬ 
larly long, of the order several years. This means that within the 
80 days of K2 observations, we would not expect to see large 
variation. With the later NITES transit however, we are starting 
to see the high amplitude TTV curve that these planets exhibit. 

The period P super depends only on A, and the period of the 
outer body. The TTV amplitude however generally shows a de¬ 
generacy between the free eccentricity Zf ree and planetary mass 
(FI2). Hadden & Fithwick|( |2014| ) break this degeneracy statis¬ 
tically, but for individual objects it can be difficult to circum¬ 
vent (although so-called synodic chopping signals can help, see 
|Nesvomy & Vokrouhlicky||2014[ )). Furthermore with our limited 
observations the current phase of the TTV curve 0 is unclear. In 
the low free eccentricity case (\Zf ree \ « |A|), 0 must be zero 
(FI2), however there is no guarantee that this is the case here. 


7.2. Transit timing variations - analysis 

We here explore the parameter space allowed by this analyti¬ 
cal representation, given the transit times we have observed. The 
analysis is most strongly constrained by the NITES transit, as 
this transit allows much more coverage of P super . The process is 
complicated by a correction term which must be added to our de¬ 
rived periods. As they are derived from a limited part of the TTV 
phase curve, they do not represent the overall ‘true’ period, in 
the sense of the mean transit interval over long timescales. This 
makes determination of P super non-trivial. Small corrections to 
the periods can change A significantly, which has a strong effect 
on the TTV period and amplitude. We circumvent this problem 
by utilising the fact that our derived periods are in fact measure¬ 
ments of the gradient of the TTV curve at the time of the K2 
observations. The sensitivity to A exhibited by the TTV period 
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Fig. 6. The x 2 surface given by the assumption that \Zf ree \ « |A|. A 
good fit is obtained for limited combinations of the two masses. Contour 
levels are shown at intervals of 0.5. 

and amplitude cancels out when calculating the gradient, allow¬ 
ing the correction to the periods to be made using only the initial 
periods. Using this, we can fit our transit times using the follow¬ 
ing process: 

1. Take values for the input parameters, Mb, M c , 7te(Zf ree ), 
/m(Zf ree ), 0 

2. Calculate the correction to make to the derived periods 

3. Use the corrected periods to find the true A and P sup er 

4. Use the input parameters and the newly corrected values to 
calculate the TTV amplitude 

5. Compare the now fully defined TTV curve to the observa¬ 
tions 

We begin with the case where \Zf ree \ « |A|, and the free 
eccentricity can be ignored. This reduces the input parameters to 
solely the two planetary masses, as 0 must also be zero in this 
case (LI2). We set to to be zero at the time of first transit - this is 
accurate to within -20 days, and given the long Psuper of several 
years, this does not need to be more accurate for this analysis. 

The parameter space of the low eccentricity case is best 
shown in Figure [b] which shows the log x 2 surface seen. The 
key points are a maximum mass for the outer planet c, of 350M®, 
found when the mass of the inner planet b goes to zero. The mass 
of planet b is less constrained (to constrain it properly using this 
analysis would require observations of planet c’s TTV curve), 
but can only take high values in the event that planet c’s mass 
drops much lower. Similarly F super , excepting very high masses 
(over ~6OOM 0 ) for planet b, is constrained to be greater than 
-480 days, and less than -3050 days. The accompanying surface 
for the ‘true’ A is shown in Figure [7J and makes clear that all 
zero eccentricity best fitting planetary mass combinations imply 
a ‘true’ A that is in fact slightly below the 3:2 resonance. As such 
while we cannot confirm this is the case from these observations 
alone, it is possible that the apparent planetary periods oscillate 
around the resonance over the course of the TTV curve. 

Extending our analysis to the case where there is significant 
free eccentricity, we can immediately constrain 0. Because the 
NITES transit arrived late rather than early, 0 must be in the 
range 0 < 0 < n. Within this range however a number of differ¬ 
ent effects can occur. We test these cases by repeating the analy¬ 
sis with the real component of Zf ree set to be 0.2, at various val¬ 
ues of 0. We hold the imaginary component at zero. Although 
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Fig. 7. As Figure [6] but showing the normalised distance to resonance 
A over the same mass ranges. In this case it can be seen that A becomes 
negative for all good fits, implying that the planets are in fact below 
the 3:2 resonance if this assumption held. Contour levels are shown at 
intervals of 0.004. 



M c (M e ) 

Fig. 8. As Figure[6] but showing the^ square surface for Re(Z free )=0.2, 
and 0 = 7r/4 over the same mass range. A degeneracy can be seen 
between positive A (top left) and negative A (down and towards the 
right). Contour levels are shown at intervals of 0.5. 

the imaginary component can affect the amplitude and phase of 
the TTV curve, as we are trialling different values of 0 it is prin¬ 
cipally the amplitude of Zf ree which is important. The results of 
these tests can be summed up simply: while M\ remains poorly 
constrained, M c can only rise above its zero eccentricity value in 
vary rare cases, and then not by much. The worst case we found 
was for 0 = 7 t/4 , where for Mb - 0 the maximum mass for M c 
was 386M 0 . The particular^ 2 square surface varies for different 
input 0 values. One example surface is shown in Figure [8j which 
also demonstrates an interesting degeneracy that arises between 
positive and negative A in that case. 

Given the poorly mapped TTV curve, the sensitivity of our 
analysis to small corrections to the planetary periods, and the 
free eccentricity degeneracy it would be premature to make de¬ 
terminations of the planetary masses at this point. We can how¬ 
ever limit both them and the corrections to the periods which 
would have to be made. As has been stated, M c < 350M 0 in the 
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zero free eccentricity case. At this point (Mb = 0, M c = 350M©) 
we also obtain the maximum period correction to make to planet 
b. This is +0.029 days to Pb in the zero free eccentricity case. 
As planet b’s mass is less constrained placing a limit on the pe¬ 
riod correction for planet c is harder, but limiting Mb to 600M© 
gives a maximum amplitude for the period correction to P c of 
-0.12 days. At the zero eccentricity case (0 = 0) we are already 
at the steepest gradient found on the TTV phase curve - as this 
sets the period correction, these amplitudes cannot go higher. 
They can however change sign (at 0 = n for example), and so 
we constrain |Correction^)| < 0.029d and |Correction(P c )| < 
0.12d immediately. The correction to Pb can be further limited 
by noting that not all of the masses which provide a fit at 0 = 0 
do so at other values for 0. In particular, when 7r/2 < 0 < n (the 
case for a negative correction to Pb), the allowed range for M c is 
much smaller (M c < 30M©), which corresponds to a maximum 
negative correction of Correction( J P^,)>-0.002. As such, the final 
limits for Pb are -0.002 < Correction^) < 0.029, where Cor¬ 
rection^) is to be made to the periods found in Table [2] This 
limits -0.011 < A < 0.013 in the extreme case, confirming that 
the system remains very close to resonance. 

At this stage it is worth stating the now better understood 
periods of these two objects. Using the period implied by our 
latest NITES transit measurement and the To from K2, we obtain 
Pb = 7.921+0003 days , and P c = 11.91 ± 0.12 days, where the 
errors are ranges rather than lcr errors, and account for TTV 
related period corrections. When predicting transit times, these 
periods and the To values of Table ^should be used. Note that 
there will also be possible TTVs of magnitude up to at least an 
hour. 


7.3. Hill stability 

Stable main sequence evolution of a two-planet system may be 
guaranteed by residing in a 3:2 MMR, although the Kirkwood 
gaps demonstrate that this MMR can instead harbour unstable 
orbits. Also, as-yet-undetected planets perturbing the 3:2 MMR 
can cause complex dynamical structures (Fuse] [2002] ) and po¬ 
tentially instability. The potential protection afforded by the 3:2 



Fig. 9. Hill stability limits for the K2-19 system. We assume M* = 
0.9 Mq, the mutual inclination between the planets is 1 degree, and a h /a c 
is within 0.1 per cent of (3 /2) 2/3 based on the planetary orbital period 
ratios. Eccentricites are measured in Jacobi coordinates. If K2-19 is Hill 
unstable, then the 3:2 MMR might act as a crucial protection mechanism 
to ensure the system’s long-term stability on the main sequence. 


of the eccentricities of both bodies (measured in Jacobi coordi¬ 
nates) does not exceed about 0.2, and (ii) that each planet is less 
massive than Jupiter. These relations help motivate the setup for 
Figure [9] The figure plots the pairs of planet masses for different 
eccentricities (measured in Jacobi coordinates) that would place 
the system on the edge of Hill stability, assuming a mutual in¬ 
clination of 1 degree and a semimajor axis ratio which is just 
0.1 per cent within (3/2) 2/3 . We generated each set of coloured 
points by sampling 600 different values of each of Mb and M c 
uniformly in log space between 10 _6 M© and 10 _2 M o , a range 
which covers both Earth masses and Jupiter masses. The plot 
axes span the entire range of masses that we sampled. 

Even if K2-19 is Hill stable, then planet c might eventu¬ 
ally escape the system or planet b might crash into the star 


MMR becomes more important when the two planets are Hill 

through Lagrange instability (Barnes & Greenberg 2006, 2007; 

unstable. Two planets are Hill stable if their orbits are guaran- 

Raymond et al. 2009; 

Kopparapu & Barnes 2010; Deck et al. 

teed to never cross; hence Neptune and Pluto are Hill unstable, 

2013 ; Veras & Mustill 

20131. No analytical Lagrange unstable 


a function of masses, semimajor axes, eccentricities and incli¬ 
nations. [Veras et al.[( |2013| outline an algorithm for computing 


the Hill stability limit; no explicit formula exists for arbitrary 
eccentricities and inclinations. 

In order for K2-19b & c to be Hill unstable, their masses 
and/or eccentricities must be sufficiently large. The mutual in¬ 
clination between the planets of just about a degree n egligibly 
affects the Hill stability limit (Veras & Armitage||2004 \ Con¬ 


straining the eccentricities and masses based on orbital periods 
alone with Hill stability is a useful exercise but requires assump¬ 
tions. A commonly-made assumption for transiting planets is 
that those planets are on circular orbits; the closer the planet is to 
the star, the better that assumption, based on tidal circularisation 
arguments. We need not make such assumptions here. 

We can use the green curves in Fig. 1 of |Veras & Ford] ( [2QT2| 
to roughly estimate Hill stability limits for K2-19. Broadly, the 
plot shows that the system will be Hill stable if both (i) the sum 


7 Just the existence of a nonzero mutual inclination between the planets 
indicates that the planets might in stead reside in an inclination-based 
6:4 MMR (see equation 1 of Milani et al. 1989]) . 


boundary is known to exist. Regardless, the 3:2 MMR may then 
provide a protection mechanism not only for Hill unstable sys¬ 
tems, but also for Lagrange unstable systems. The upcoming 
space mission PLATO ( [Rauer et al.|2014| ) will provide accurate 
enough stellar age constraints to potentially detect a decreasing 
trend in stable multi-planet systems with time due to Lagrange 
instability ( Veras et al .~|2015| ). 


8. Conclusion 

We have presented and validated K2-19b & c, a system of two 
Neptune sized planets orbiting close to the 3:2 MMR, via tran¬ 
sit observations with K2 and the NITES telescope and analysis 
using PASTIS. The inner, larger planet b shows high amplitude 
TTVs, and will likely show larger amplitudes when further tran¬ 
sits are observed. The outer, smaller planet c can be expected to 
show even greater TTVs (scaled up by the mass ratio of the two 
planets), although these have not yet been observed. The precise 
ephemeris of these planets is still in doubt, see Sections]?] and [ 7 ] 
for the fit values and limits to the possible corrections to them in 
the near but not in resonance case. 
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Future observations of K2-19 have the potential to lead to 
interesting discoveries. The system is bright enough to observe 
from the ground, leading to great potential for future work. The 
observation of more transits of either planet will lead to fully 
characterised TTV phase curves, as well as possibly being able 
to fully solve for the planetary masses via full dynamical anal¬ 
ysis. Radial velocity observations have the potential to indepen¬ 
dently characterise the planetary masses, providing a window on 
the discrepancies that sometimes exist between radial velocity 
and TTV derived masses. We expect that this system is one of 
many more which will arise from the K2 mission. 
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Appendix A: PASTIS supplements 



Fig. A.l. Sketchs of the different scenarios considered for the PASTIS validation of the two planets in the K2-19 system. The sizes are not 
proportional to their real values but note the relative size of the orbits between the 8-d and 12-d signals. Note also the orbit of the contaminant 
system in scenarios S\ to S 5 , which marks the difference between these and the unbound scenarios <S 6 to <$io- 
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Table A.l. Priors used in the PASTIS analyses: Ilia, b ) represents a Uniform prior between a and b\ Nip , cr 2 ) represents a Normal distribution with a mean of p and a width of cr 2 ; !P(cr; 
represents a Power Law distribution with an exponent a computed between x min and x max ; Vo; x min ; x max ) represents a double Power Law distribution with an exponent a\ computed 

between x min and x 0 and an exponent a 2 computed between x 0 and x max ; and finally 5(a, b) represents a Sine distribution between a and b. 


Parameter 5o 5i S 2 S 3 S 4 S 3 S 6 57 5s Sg 5io 

Target star parameters 

Effective temperature T eff [K] .A/ r (5230; 417) . 

Surface gravity log g [cm.s -2 ] .7V(4.47; 0.83) . 

Iron abondance [Fe/H] [dex] .Af(0.38; 0.23) . 

Distance d t [pc] . !P(2.0; 10; 10000) . 


Primary blend star parameters 

Stellar mass M* : [Mq] 

Iron abondance [Fe/H] +1 [dex] 
Ager*, [Gyr] 

Distance d+ l [pc] 

n/a 

n/a 

n/a 

n/a 

n/a 

n/a 

n/a 

n/a 

n/a 

n/a 

n/a 

n/a 

n/a 

n/a 

n/a 

n/a 

.. ?M-1.3;-2.3;0.5;0.1;2) .. 

n/a . 

n/a . 

n/a . 

.<Z/(-2.5;0.5) . 

.<1/(0; 13.8) . 

. P(2.0; 10; 10000) . 


Secondary blend star parameters 

Stellar mass M* 2 [Mq] 

n/a 

n/a 

n/a 

n/a 


. ^(-1.3;-2.3; 0.5; 0.1; 2) . 

. n/a n/a 

n/a 

8 -d planet parameters 

Planetary radius [R j up \ 



.. <1/(0; 2.2) .. 



n/a <1/(0; 2.2) 

n/a .<1/(0; 2.2) .. 



12 -d planet parameters 

Planetary radius R c p [R >p ] .^(0; 2.2). n/a ^/(0; 2.2) n/a .^(0; 2.2) 


8 -d orbit parameters 


Orbital period P b [d] 

Epoch T b [BJD - 2450000] 
Orbital inclination i b [°] 


A/ r (7.919454; 8 x 10 -3 ) 

. AA(6813.38345;0.03) . 
.5(70; 90). 


12 -d orbit parameters 


Orbital period P c [d] . N(l 1.90701; 4 x 10“ 2 ) 

Epoch T c 0 [BJD - 2450000] .7V(6817.2759; 0.1) . 

Orbital inclination i c [°] .5(70; 90) . 
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Table A.l. Continued. 

Parameter So S\ S 2 S 3 S 4 S 5 S6 S 7 Sg S\o 

Kepler data parameters 

Contamination .W(0; 0.02). 

Out-of-transit flux . < Z/(0.9; 1.1) . 

Jitter .*Z/(0; 0.1) . 

Spectral Energy Distribution parameter 

Jitter [mag] .^/(0; 1). 
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Table A.2. Posterior distributions of the PASTIS analyses. 


Parameter 


<Si 

s 2 

S 3 

s 4 

<$5 

^6 

<$7 

S 8 

S 9 

<$10 

Target star parameters 

Effective temperature T ef f [K] 
Surface gravity logg [cm.s -2 ] 

Iron abondance [Fe/H] [dex] 
Distance d t [pc] 

5562±41 

4.509 

+0.026 

-0.051 

0.07+0.16 

329!“ 

5774+51 

4.397 

±0.036 

0.24+0.17 

501+18 

5756+53 

4.426 

+0.033 

-0.052 

0.03+0.18 

470+19 

5732+58 

4.358 

+0.057 

0.03+0.20 

499+21 

5550+130 

4.527 

+0.037 

-0.061 

-0.08+0.15 

452+23 

5752+55 

4.332 

+0.033 

-0.050 

0.14+0.17 

521+24 

5487+61 

4.508 

+0.059 

0.10+0.19 

329+25 

5308+98 

4.27 

+0.20 

0 33 +(m 
u '~’~’-0.18 

443 + 96 

5490+66 

4.508 

+0.058 

0.08+0.19 

329+27 

5333+210 

4.441 

+0.072 

-0.190 

0.29+0.18 

730 +480 

5368+84 

2.88 

+0.35 

0-3 8!^ 
3800+1900 

Primary blend star parameters 












Stellar mass [Mq] 

n/a 

0.914 

0.862 

0.858 

0.265 

0.883 

1.16 

0.933 

1.00 

0.824 

0.90 



+0.048 

+0.059 

+0.058 

+0.026 

+0.048 

+0.18 

+0.110 

-0.068 

+0.26 

-0.15 

+0.079 

+0.11 

Iron abondance [Fe/H]*j [dex] 

n/a 

n/a 

n/a 

n/a 

n/a 

n/a 

-1.8+0.5 

-1.5+0.7 

o o 
+ l 

OO 

-0.6-; 3 

-1.8+0.6 

Age t*j [Gyr] 

n/a 

n/a 

n/a 

n/a 

n/a 

n/a 

1 19 +16 

-0.76 

U4 

to 

o 4^ 

1 5 +3 - 8 
A * J -1.1 

5.6+4.4 

9 1 +4.4 

A -1.6 

Distance d+ ] [pc] 

n/a 

n/a 

n/a 

n/a 

n/a 

n/a 

3200 

820 

2260 

336 

880 








+1300 

+320 

-120 

+1300 

-720 

+ 190 
-39 

+240 

-150 


Secondary blend star parameters 

Stellar mass M* 2 [M©] n/a n/a n/a n/a 0.928 0.103 0.186 0.112 n/a n/a n/a 

+0 056 +0005 +0 064 +0016 

=nvj.u~>vj -0.002 -0.009 


8 -d planet parameters 












Planetary radius [Rj up ] 

0.661 

+0.053 

-0.033 

0.960 

+0.040 

1.087 

+0.053 

1.141 

+0.050 

0.906 

+0.037 

n/a 

0.677 

+0.060 

n/a 

0.678 

+0.063 

0.661 

+0.26 

-0.069 

1.268 

+0.12 

12 -d planet parameters 












Planetary radius R^ [Rj up ] 

0.385 

+0.031 

-0.018 

0.679 

+0.029 

0.521 

+0.025 

-0.016 

0.669 

+0.028 

n/a 

0.580 

+0.037 

-0.021 

n/a 

0.60 

+0.17 

1.65 

+0.32 

0.99 

+0.55 

-0.33 

0.743 

+0.070 

8 -d orbit parameters 

Orbital period P b [d] 

Epoch T b [BJD - 2456000] 

Orbital inclination i b [°] 

7.919491 

±8.3xlO -5 

813.38351 

±3.8xl0 -4 

88.7!“ 

7.919510 

+9.0xl0 -5 

813.38347 

+4.1X10" 4 

87.8+0.2 

7.919527 

+9.0xl0 -5 

813.38340 

+3.9xl0 -4 

89.7+0.3 

7.919525 

+9.1xl0 -5 

813.38341 

+4.0xl0 -4 

89.1+0.3 

7.919502 

+8.6xl0 -5 

813.38348 

+4.1xl0 -4 

88.9+0.6 

7.919520 

+9.0xl0 -5 

813.38345 

+4.2xl0 -4 

89.5+0.4 

7.919494 

+8.3xl0 -5 

813.38350 

+3.8xl0 -4 

88.7+0.5 

7.919486 

+8.6xl0 -5 

813.38358 

+4.0xl0 -4 

89.4+0.5 

7.919486 

+8.7xl0 -5 

813.38351 

+3.8xl0 -4 

88.7+0.5 

7.919477 

+8.7xl0 -5 

813.38355 

+3.9xl0 -4 

89.4+0.5 

7.919501 

+9.0xl0 -5 

813.38351 

+4.1xl0 -4 

89.1+0.2 


Armstrong et al.: K2-19b & c 












Article number, page 18 of|l8| 


Table A.2. Continued. 


Parameter 

•So 

Si 

<s 2 

<S 3 

S 4 

<$5 

^6 

<s? 

<s 8 


<$10 

12 -d orbit parameters 

Orbital period P c [d] 

Epoch T c 0 [BJD - 2456000] 

Orbital inclination i c [°] 

11.90664 
+4.9xl0 -4 
817.2757 
+ 1.4xl0 -3 
89.4+0.5 

11.90634 
±4.9x10 -4 
817.2758 
±1.4xl0 -3 
89.81“ 

11.90661 

±4.7xl0 -4 

817.2758 

±1.4xl0“ 3 

88.71“ 

11.90642 
+5.0xl0 -4 
817.2758 
+ 1.4xl0 -3 
89.8+0.2 

11.90636 

±5.1xl0 -4 

817.2760 

±1.5xl0 -3 

89.8±0.2 

11.90730 
+5.0xl0 -4 
817.2735 
+ 1.5xl0 -3 
88.1+0.2 

11.90650 

±5.2xl0 -4 

817.2760 

±1.5xlO -3 

89.4±0.6 

11.90721 
+4.9xl0 -4 
817.2739 
+ 1.5xl0 -3 
87.9+1.2 

11.90650 

±4.9xl0 -4 

817.2760 

±1.5xl0 -3 

89.4±0.6 

11.90659 

±5.1xl0 -4 

817.2760 

±1.4xl0 -3 

88.912;? 

11.90637 
+5.0xl0 -4 
817.2758 
+ 1.5xl0 -3 
89.7+0.3 

Kepler data parameters 

Contamination 

Out-of-transit flux 

Jitter [ppm] 

-0.001 

±0.021 

1.0000223 

+9.0xl0 -6 

59+23 

-0.003 

+0.020 

1.0000228 

+9.2xl0 -6 

63+21 

-0.003 

+0.020 

1.0000236 

+9.6xl0 -6 

67+20 

-0.004 

+0.021 

1.0000238 

+9.4xl0 -6 

68+19 

-0.002 

+0.020 

1.0000101 

+9.2xl0 -6 

68+20 

0.001 

+0.021 

1.0000320 

+9.3xl0 -6 

68+20 

-0.002 

+0.021 

1.0000228 

+9.3xl0 -6 

60+23 

-0.001 

+0.021 

1.0000315 

+9.9xl0 -6 

60+22 

-0.001 

+0.020 

1.0000234 

+9.4xl0 -6 

60+23 

-0.001 

+0.020 

1.0000221 

+8.9xl0 -6 

58+24 

-0.002 

+0.020 

1.0000230 

+9.5xl0 -6 

65+22 

Spectral Energy Distribution parameter 











Jitter [mag] 

0.04!°f 

0.04!°f 

0.0412” 

0.0412” 

0.0412” 

0.0412;” 

0.0412;” 

0.0412” 

0.0412” 

0.0412” 

0.0412;” 
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